Pattern formation of microtubules and motors: inelastic interaction of polar rods 
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We derive a model describing spatio-temporal organization of an array of microtubules interacting 
via molecular motors. Starting from a stochastic model of inelastic polar rods with a generic 
anisotropic interaction kernel we obtain a set of equations for the local rods concentration and 
orientation. At large enough mean density of rods and concentration of motors, the model describes 
orientational instability. We demonstrate that the orientational instability leads to the formation 
of vortices and (for large density and/or kernel anisotropy) asters seen in recent experiments. 

PACS numbers: 87.16.-b, 05.65.+b,47.55.+r 



One of the most important functions of molecular mo- 
tors (MM) is to organize a network of long filaments (mi- 
crotubules, MT) during cell division to form cytoscele- 
tons of daughter cells hi. A number in vitro experiments 
were performed 0, EI H IE El to study interaction of 
MM and MT in isolation from other biophysical processes 
simultaneously occurring in vivo. At large enough con- 
centration of MM and MT, the latter organize in asters 
and vortices depending on the type and concentration of 
MM. 

After MM binds to a microtubule at a random posi- 
tion, it marches along it in a fixed direction until it un- 
binds without appreciable displacement of MT. If a MM 
binds to two MTs, it can change their mutual position 
and orientation significantly. In Ref. [l| , the interaction 
of rod-like filaments via motor binding and motion has 
been studied, and patterns resembling experimental ones 
were observed. In a phenomenological model for the 
MM density and the MT orientation has been proposed. 
Ref. generalized this model by including separate den- 
sities of free and bound MM, as well as the density of MT. 
They found the transition from asters to vortices as the 
density of MM is increased, in disagreement with exper- 
imental evidence that the asters give way to vortices 
with decreasing the MM concentration. A phenomenolog- 
ical flux- force relation for active gels was applied in |10| . 
While vortex and aster solutions were obtained, an anal- 
ysis of that model is difficult because of a large number 
of unknown parameters. In Ref. a set of equations 
for MT density and orientation was derived by averaging 
conservation laws for MT probability distribution func- 
tion. However, this model does not exhibit orientation 
transition for the homogeneous MT distributions. 

Here we derive a model for the collective spatio- 
temporal dynamics of MTs starting with a master equa- 
tion for interacting inelastic polar rods. Our model differs 
from the transport equations 01 m that it maintains the 
detailed balance of rods with a certain orientation. The 
model exhibits an onset of orientational order for large 
enough density of MT and MM, formation of vortices and 
then asters with increase in the MM concentration. 



MMs enter the model implicitly by specifying the in- 
teraction rules between two rods. Since the diffusion of 
MMs is about 100 times higher than that of MTs, at the 
first step we neglect spatial variations of the MM density. 
While variable MM concentration affects certain quanti- 
tative aspects ,7], our analysis captures salient features 
of the phenomena. All rods are assumed to be of length 
I and diameter d -C I, and are characterized by the its 
center of mass r and orientation angle 4>. 

Maxwell model. Consider the orientational dynamics 
only and ignore the spatial coordinates of interacting 
rods (an analog of the Maxwell model of binary colli- 
sions in kinetic theory of gases, see e.g. ^3)- Since 
motor residence time on MT (about 10 sees) is much 
smaller than the characteristic time of pattern formation 
(about 1 hour), we model MM-MT inelastic interaction 
by an instantaneous collision in which two rods change 
their orientations: 
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after the col- 



where cf>\ 2 are orientations before and 
lision, and 7 characterizes inelasticity of collisions, and 
102 — 4>\ I < 0o < t. The angle between two rods is 
reduced after the collision by a factor 2j — 1. 7 = 
corresponds to a totally elastic collision (the rods ex- 
change their angles, this case does not have counter- 
part in MM-MT interaction context) and 7 = 1/2 corre- 
sponds to a totally inelastic collision: rods acquire iden- 
tical orientation 4>i 2 = (01 + ^2)/^ ( see Fig. H a )- Here 
we assume that two rods only interact if the angle be- 
tween them is less than <pQ. Because of 27r-periodicity, 
we have to add the rule of collision between two rods 
with 2ir — <fio < |<^2 — (f>i\ < 2tt. In this case we have to 
replace <fi{ a — > <fi{ a + 7r, 2 ' a — * 2 ' a — n in Eq. (|TJ. In the 
following we will only consider the case of totally inelastic 
rods (7 = 1/2) and <po = the generalization for arbi- 
trary 7 and <))q is straightforward [lflj . The probability 
P((j>) obeys the following master equation 
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d t P(fa = D r 8%P(<l>) + g 



)P{fo) (2) 



x [S(4> - fa/2 - (hi 2) - 5^ - fa)} + 9 / dfadfa 

Jc 2 

xP(fa)P(fa)[5{<t> ~ <h/* ~ - tt) - 5{<i> - fa)} 

where g is the "collision rate" proportional to the num- 
ber of MM, the diffusion term oc D r describes the ther- 
mal fluctuations of rod orientation, and the integration 
domains C\ , C 2 are shown in Fig(IJi. Changing variables 
t — ► D r t, P — > gP/D r , w = fa — fa , one obtains 



d t P(fa = dlP(fa 



dw 



x [P((j> + w/2)P(4> - w/2) - P{faP{(j) - w)\ (3) 

The rescaled number density p = P{fa t)d(f> now is 
proportional to density of rods multiplied by the density 
of motors. Let us consider the Fourier harmonics: 



Pk = (e 
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The zeroth harmonic Pq = p/2n = const, and the real 
and imaginary parts of Pi represent the components t x = 
(cos fa ,T y = (sin fa of the average orientation vector t, 
t x + iT y — P*. Substituting (@J into Eq.® yields: 

Pk + {k 2 + p)P k = 2tt Pk- m P m S[Tvk/2 - rrm] (5) 



(here S(x) = sinx/x). Due to the angular diffusion term, 
the magnitudes of harmonics decay exponentially with 
|fc|. Assuming P^ = for \k\ > 2 one obtains from Eq.® 
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P 1 =P P 1 2(4-tt)--P 2 P* 



P 2 + 4P 2 = -P P 2 2ir + 2-kP{ 



(6) 
(7) 



Since near the instability threshold the decay rate of 
P 2 is much larger than the growth rate of P\, we can 
neglect the time derivative P 2 and obtain P 2 = AP 2 with 
A = 2ir(p + 4) _1 and arrive at: 



et-Aq\t\ 2 t 



(8) 



withe = p(4 7 r- 1 -l)-f 0.273p-I and A = 8A/3. For 
large enough p > p c = 7r/(4 — n) w 3.662, an ordering 
instability leads to spontaneous rods alignment. This 
instability saturates at the value determined by p. Close 
to the threshold Aq m 2.18. Fig. [3 shows stationary 
solutions P(fa) obtained from Eq. ©. As seen from the 
Inset, the corresponding values of |r| are consistent with 
the truncated model (JSJ up to p < 5.5. 

To describe the spatial localization of interactions we 
introduce the probability distribution P(r, fa t) to find 
a rod with orientation <j) at location r at time t. The 
master equation for P(r, fa t) can be written as 



d t P(r, fa = dlP(v, fa + diDijdjPir, , 



xP(r 1 ,4> + w/2)P(r 2 ,(j)~w/2)5 [ Tl+Y2 _ r 



/4>o 
dw [W(ri, r 2 , <p + w/2, <f> - w/2) 
-00 

-W(r u r 2 , fa<f>- w)P(r 2 , 0)P(r x , <f> - w)S (r 2 - r)] (9) 



where we performed the same rescaling as in Eq.® 
and dropped argument t for brevity. The first two 
terms in the r.h.s. of describe angular and 
translational diffusion of rods with the diffusion ten- 
sor = jj- [D^riinj + Dx(Sij — niUj)) . Here n = 
(cos(fa), sm(fa)). D r ,D\\,D± are known in polymer 
physics: Dy = *f^,D± = ^,D r = where 
£i|,£_L,£r are corresponding drag coefficients. For rod- 
like molecules, £|| = 2Trn s l/\og(l/d); £_l = 2£||;£ r « 
nr) a l 3 /31og(l/d) where n s is shear viscosity Q. 

The last term of Eq. © describes MM-mediated inter- 
action of rods. We assume that after the interaction, the 
two rods acquire the same orientation and the same spa- 
tial location in the middle of their original locations. The 



interaction kernel W is localized in space, but in general 
does not have to be isotropic. On the symmetry grounds 
we assume the following form (we assume 2D geometry 
and neglect higher-order anisotropic corrections): 



W = -rj— exp 

Z TT 



(ri - r 2 ) 2 
b 2 



(1 +/?(!•! -r 2 ) • (m -n 2 )) 



with b m I = const. This form implies that only 
nearby MTs interact effectively due to MMs. The 0(0) 
anisotropic term describes the dependence of the cou- 
pling strength on the MT mutual orientation: "diverg- 
ing" polar rods (such as shown in Fig. Q] a ) inter- 
act stronger than "converging" ones. This is the sim- 
plest term yielding non-trivial coupling between density 
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and orientation. We perform Fourier expansion in 4> 
and truncate the series at |n| > 2, 2-kP gives the lo- 
cal number density p(r,t), and P±\ the local orienta- 



tion r(r, t). Omitting calculations (see |l3J), rescaling 
space by Z, and introducing dimensionless parameters 
B = b/L H = /31B 2 , we arrive at 



d t p = 
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d t r = 
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[3V • (rV 2 p - pVV) + 25, (c^r 4 - bpdjTj)] 
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The last two terms in Eqs. . (|1 1|1 are linearized near 
the mean density po = (p). The last term in Eq. I|10(l 
regularizes the short-wave instability when the diffusion 
term changes sign for p > pb = 1/4_B 2 . This instability 
leads to strong density variations associated with forma- 
tion of MT bundles, see Fig. 

Aster and vortex solutions. If B 2 H <C 1, the density 
modulations are rather small, and Eq. Ijlljl for orien- 
tation r decouples from Eq. (|10[1 . It is convenient to 
rewrite Eq. (|ll|l for complex variable ip — T x +ir y in polar 
coordinates r, 9: i\) = F(r) exp[i8 + i<p(r)] where the am- 
plitude F(r) and the phase <p(r) are real functions. For 
the aster solution (p(r) — and for the vortex <p{r) =/= 0. 
Asters and vortices can be examined in the framework of 
one-dimensional problem for V = \/AoF(r) exp[iip(r)]: 

d t V = DxArV + D 2 A r V* + (l - \V\ 2 ) V 

-H ^aiFReV r F + a 2 d r VReV + a2 ^ mV ^ ( 12 ) 

where A r = d 2 + r~ 1 d r — r~ 2 , V r = d r + r -1 , D\ = 
1/32 + /9 S 2 /47r, D 2 = 1/192, ai = (tt - 8/3)/v^ « 
0.321, a2 = 8/3yj4o ~ 1.81, and we rescaled time t — > t/e 
and space by r — > r/^/e. The aster and vortex solutions 
for certain parameter values obtained by numerical in- 
tegration of Eq. (|12|) are shown in Fig. 0] Vortices 
are observed only for small values of H and give way 
to asters for larger H. For H = 0, Ea. (|12|) reduces to 
a form which was studied in [l5| . It was shown in [l5| 
that the term A r V* favors vortex solution (<p = tt/2). 
In contrast, the terms proportional to H select asters. 
Increasing H leads to gradual reduction of ip, and at a 
finite -ffo(po) 4'{ r ) = 0? i- e - l ne transition from vortices 
to asters occurs. For < H < Hq, the vortex solution 
has a non-trivial structure. As seen in Fig. 0] the phase 
ip — > for r oo, i.e. vortices and asters become indis- 
tinguishable far away from the core. 

The phase diagram is shown in Fig. The solid line 
Ho(po) separating vortices from asters is obtained from 
solution of the linearized Eq. Ijl2(l by tracking the most 
unstable eigenvalue A of the aster. For this purpose the 
solution to Eq. i|12|) wa s sought in the form V = F + 



iwexp(At), where real w obeys L = Xw with operator 

L = DA r + (1 - F 2 - ai HV r F) - a 2 HFV r (13) 

(D = D\ — D 2 ). Eq. 113|) was solved by the matching- 
shooting method. The dashed line corresponds to the 
orientation transition limit po — p c . The lines meet in 
a critical point H c = Ho(p c ) above which vortices are 
unstable for arbitrary small e > 0. The phase diagram is 
consistent with experiments, see Ref. [fj: for low value 
of kernel anisotropy H < H c (possibly corresponding to 
kinesin motors) increase of the density po first leads to 
formation of vortices and then asters. For H > H c (pos- 
sibly corresponding to Ned) only asters are observed. 

For H 7^ well-separated vortices and asters exhibit 
exponentially weak interaction. For asters it follows from 
the fact that L is not a self-adjoint operator. Null-space 
of U exponentially decays at large r, w ~ exp[— r/Lo] 
with the screening length in the original length units 
Lo = D/a 2 H^/e (see Q). 

We studied the full system (|10fl . (|ll|) numerically. In- 
tegration was performed in a two-dimensional square 
domain with periodic boundary conditions by quasi- 
spectral method. For small H we observed vortices and 
for larger H asters, in agreement with the above analy- 
sis. As seen in Fig. \5\ asters have unique orientation of 
the microtubules (here, towards the center). Asters with 
opposite orientation of r are unstable. In large domains 
asters form a disordered network of cells with a cell size 
of the order of Lq. Neighboring cells are separated by 
the "shock lines" containing saddle-type defects. Start- 
ing from a random initial condition we observed initial 
merging and annihilation of asters. Eventually, annihi- 
lation slows down due to exponential weakening of the 
interaction of asters. 

We derived continuous equations for the evolution of 
MT concentration and orientation. We found that ini- 
tially disordered system exhibits an ordering instability 
similar to a nematic phase transition in ordinary poly- 
mers at high density. The important difference is that 
here the ordering instability is mediated by MM and can 
occur at arbitrary low densities of MT. At the nonlin- 
ear stage, the instability leads to experimentally observed 
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FIG. 2: Stationary solutions P((j>) for different p. Inset: the 
stationary value of t vs p obtained from the Maxwell model 
©, dashed line - truncated model 



FIG. 1: a - sketch of motor-mediated two-rod interaction for 
7 = 1/2, b - integration regions Ci,2 for Eq.J5J. 
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FIG. 3: Phase boundaries obtained form the linear stability 
analysis of aster solution for B 2 = 0.05, dashed line shows 
bundling instability limit po > pb — 5. Inset: Position of 
critical point H c vs B at po — 4.5. 
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FIG. 4: Stationary vortex and aster solutions t x + it v — 
F(r) exp{i8 + iip(r)] to Eq. JEJ, for po = 4, B 2 = 0.05. 




FIG. 5: Orientation r for vortices (H — 0.006, left) and asters 
(H = 0.125, right) obtained from Eqs. pOllip . Color code 
indicates the intensity of |t| (red corresponds to maximum 
and blue to zero), B 2 = 0.05, po = 4, domain of integration 
80 x 80 units, time of integration 1000 units. 



